function [Xi,Ti,Ui,cell_type] = el_resample(grid,base,U,varargin)
% [Xi,Ti,Ui] = el_resample(grid,base,U)
% [Xi,Ti,Ui] = el_resample(grid,base,U,k)
% [Xi,Ti,Ui] = el_resample(grid,base,U,k,l)
% [Xi,Ti,Ui] = el_resample(grid,base,U,k,l,type)
%
% Resample the solution U (stored as two dimensional array) on a local
% mesh for each element. The two optional arguments k and l are the
% global and the local element degrees. If both are absent, base.k is
% used; if l is absent, l=k is used.
%
% If type=['hybrid'|'flux] is also present, it is assumed that the
% solution U is defined on the sides of the grid, or as a vector
% variable, and base.e_s or base.o_s are used as basis.
%
% Note: the local mesh is defined in "ResampleMesh.octxt"

 if(nargin<4)
   k = base.k;
 else
   k = varargin{1};
 end
 if(nargin<5)
   l = k;
 else
   l = varargin{2};
 end
 hyb = 0; flx = 0;
 if(nargin>=6)
   hyb = strcmp(varargin{3},'hybrid');
   flx = strcmp(varargin{3},'flux');
 end
   
 locmeshes = load("ResampleMesh.octxt","meshes");
 if(k>size(locmeshes.meshes,2))
   error(['Global degree ',num2str(k),' not supported.']);
 endif
 if(l>size(locmeshes.meshes,3))
   error(['Local degree ',num2str(l),' not supported.']);
 endif

 locmesh = locmeshes.meshes(grid.d-hyb,k,l);
 p = locmesh.p;
 t = locmesh.t;

 % 1) evaluate the basis functions on the local mesh
 np = size(p,2);
 nt = size(t,2);
 if(hyb)
   PHI = zeros(base.nk,np);
   for i=1:base.nk
     PHI(i,:) = ev_pol(base.e_s{i},p);
   end
 elseif(flx)
   PHI = zeros(grid.d,base.mk,np);
   for i=1:base.mk
     for j=1:grid.d
       PHI(j,i,:) = ev_pol(base.o_s{j,i},p);
     end
   end
 else
   PHI = zeros(base.pk,np);
   for i=1:base.pk
     PHI(i,:) = ev_pol(base.p_s{i},p);
   end
 end

 % 2) resample all the elements (sides if hyb)
 if(hyb)
   no = grid.ns; % number of "objects"
 else
   no = grid.ne;
 end
 Xi = zeros(grid.d,np,no);
 Ti = t; % the connectivity is always the same
 if(flx)
   Ui  = zeros(grid.d,np,no);
   Ui1 = zeros(       np,no);
 else
   Ui = zeros(       np,no);
 end
 idx = ones(1,np);
 if(hyb)
   o = grid.s;
 else
   o = grid.e;
 end
 for io=1:no
   Xi(:,:,io) = o(io).b * p + o(io).x0(:,idx);
 end
 if(flx)
   for i=1:grid.d
     PHI1 = permute(PHI(i,:,:),[2,3,1]);
     for io=1:no
       Ui1(:,io) = PHI1' * U(:,io);
     end
     Ui(i,:,:) = permute(Ui1,[3,1,2]);
   end
   % include the Piola transformation
   for io=1:no
     piola = o(io).b/o(io).det_b;
     Ui(:,:,io) = piola*Ui(:,:,io);
   end
 else
   for io=1:no
     Ui(:,io)   = PHI' * U(:,io);
   end
 end

 % 3) Paraview specific information
 if(nargout>3)
   cell_type = locmesh.PV_cell_type;
 endif

return

